TOPICAL REVIEW 



A stochastic approach to open quantum system 

R. Bielei'2 and R. D'Agosta^ ^ 

^Institut fiir Theoretische Physik, Technische Universitat Dresden, D-01062 
Dresden, Germany 

^ ETSF Scientific Development Center, Departamento de Fi'sica de Materiales, 

Universidad del Pais Vasco, E-20018 San Sebastian, Spain 

3 IKERBASQUE, Basque Foundation for Science, E-48011, Bilbao, Spain 

E-mail: treborB@gmx.de, roberto.dagosta@ehu.es 



Abstract. Stochastic methods are ubiquitous to a variety of fields, ranging 
from Physics to Economy and Mathematics. In many cases, in the investigation 
of natural processes, stochasticity arises every time one considers the dynamics 
of a system in contact with a somehow bigger system, an environment, that is 
considered in thermal equilibrium. Any small fluctuation of the environment has 
some random effect on the system. In Physics, stochastic methods have been 
applied to the investigation of phase transitions, thermal and electrical noise, 
thermal relaxation, quantum information, Brownian motion etc. 

In this review, we will focus on the so-called stochastic Schrodinger 
equation. This is useful as a starting point to investigate the dynamics of open 
quantum systems capable of exchanging energy and momentum with an external 
environment. We discuss in some details the general derivation of a stochastic 
Schrodinger equation and some of its recent applications to spin thermal transport, 
thermal relaxation, and Bose-Einstein condensation. We thoroughly discuss the 
advantages of this formalism with respect to the more common approach in terms 
of the reduced density matrix. The applications discussed here constitute only a 
few examples of a much wider range of applicability. 
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1. Introduction 

In the scientific investigation of a natural system, usually the first step of the modelling 
process is to consider the system as closed and isolated. The system then evolves 
following its natural laws. The direct observation of the system -our experience- helps 
us in identifying the natural laws that many systems obey. An abstraction process 
allows us to formulate a coherent set of principles we believe are valid for a wide 
class of physical systems, we therefore create a theory [T]. However, from the same 
experience we know that no system is isolated or closed. For as small as we can think 
of, there is always some leaking or coupling, let alone the observation process, that does 
not allow for a complete decoupling of the dynamics of the system from the external 
environment. One of the major challenges of theoretical physical modelling has been 
to investigate this coupling and its effects on the "natural" system. This situation 
is even more striking for a quantum mechanical system. Here, the interaction of the 
quantum system with an external device can have dramatic effects that are neither 
adiabatic nor small. A text-book example is the measurement process: we might 
think of the coupling between the quantum mechanical system and the measurement 
tool as small as possible, but the effects on the dynamics can be as large as changing 
completely the state of the quantum system. This issue is so fundamental that is 
in fact one of the distinguishing factors between a classical and quantum mechanical 
theory. 

If this is the situation, the investigation of real quantum systems might appear 
hopeless: the coupling with an external environment will destroy any information we 
have on the system itself replacing it with some kind of forced dynamics. However, 
we have to realize that in many occasions the coupling between the system and -for 
example- a thermal bath has some random nature. From this point of view, it appears 
natural in order to describe the dynamics of our quantum mechanical system coupled 
to an external environment, to use a stochastic approach that will give some "average 
behaviour." An example from classical mechanics is the Brownian motion where a 
large particle in a liquid is hit multiple times by the smaller particles forming the 
liquid. For this reason, the large particle appears to be "suspended" in the liquid. 
At a first observation it seems to be stationary, but a closer look reveals that it is 
in constant random motion. The investigation of this motion brought Langevin to 
search for an equation of motion for the large particles following the result obtained 
by Einstein [2 : Langevin modelled the force from the liquid particles as a small 
stochastic force with zero mean and with certain correlation properties, related to 
the temperature of the liquid [3]. The success of the theory essentially opened the 
field of stochastic equations. Another interesting aspect of the Langevin work must 
be pointed out. Normally, the force that determines the dynamics in the Newton 
equation depends solely on the position of all the particles forming the system. If 
this is the case, we usually say that the system is "closed" and the dynamics is solely 
determined by the initial conditions. There is a general consensus that, enlarging the 
system under investigation, one in principle could "close" the equation of motion and 
obtain a force that is only a function of the position of all the particles in the system 
and the environment. It is clear, however, that in many cases, this is not what we 
want to do. On the one hand, the number of degrees of freedom is so large to make any 
attempt at solving the problem beyond hope, while on the other hand we deal with a 
large amount of information that is essentially useless for describing the dynamics we 
are interested in. Again we can revert to the classical Brownian motion for a simple 
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example: to close the equation of motion we need to include the dynamics of all the 
particles in the liquid, surely a large number of extra degrees of freedom. Langevin 
was able to "fold" this extra degrees of freedom into the stochastic term in his effective 
equation of motion. Further development have established a few "standard" ways to 
perform this folding and nonetheless try to maintain the important features of the 
liquid dynamics. 

In describing the dynamics of a quantum mechanical system coupled to an 
external environment we will follow a similar approach. Our starting point will 
be the Schrodinger equation for the system and the environment. We will then 
operate a selection of the relevant degrees of freedom -obviously a utterly arbitrary 
step especially if the system and the environment are made of the same interacting 
particles- and integrate out those that are deemed irrelevant. In this way, we obtain 
an effective Schrodinger equation, that will be stochastic in nature, for the "state" of 
the system. The theory we will be developing in this way is similar to the formalism 
for the density matrix. There one usually, starting from the von Neumann equation for 
the total density matrix, derives an equation of motion for a reduced density matrix, 
a master equation, that entails the physical information of the subsystem of interest. 
We will show how it is possible to derive such an equation of motion for the reduced 
density matrix from the stochastic Schrodinger equation. For this reason, the latter 
has been seen as the unraveling of a master equation. 

In this review we will concentrate on stochastic Schrodinger equations that 
simulate the average behavior of a variety of condensed matter systems interacting 
with their environments. In doing so, we will build an ensemble of states and to obtain 
any physical quantity we will average over this ensemble. For this reason, we do not 
require that a single "state" of the system describes a physical quantum trajectory: for 
example, the state might be not normalized at each instance of time, and moreover the 
time evolution might change the normalization. On the other hand, in the quantum 
theory of measurement, the case where the bath is under continuous observation with 
some type of measurement device has been investigated [H [5]. This analysis leads 
to stochastic Schrodinger equations whose solutions are single system trajectories, so- 
called "true" quantum trajectories [S]. The theory is of great importance for designing 
feedback control on open quantum systems [3] , to exploit the localization property to 
reduce the number of basis states needed to represent the state vector [7] , or to monitor 
the state of a Bose-Einstein condensate IS [9]. However, whether the individual paths 
of the stochastic Schrodinger equation are true or not does not affect the validity of 
the average results we are interested in. The interested reader could possibly start 
from [S] to explore the development of the continuous monitoring theory. 

One of the questions we will try to address is what are the physical conditions for 
the establishment of a steady state, and if this steady state corresponds to any know 
thermal equilibrium between the system and the environment. We will focus on the 
case in which there is no particle exchange between the system and the environment, 
the latter then representing a thermal bath able to supply energy and momentum 
to the system. For this reason we will talk about thermal relaxation of the system 
towards some equilibrium. There are in the literature a few examples of equations 
of motion for the state of the system that are build to describe the relaxation of the 
system towards a steady state [lUl [H] . While some of them have been widely used 
to investigate the relaxation towards a steady state, they are usually not derived but 
rather assumed due to some sought characteristics of the dynamics they impose. We 
refer the interested reader to the available literature [TUl HI [12] for a more complete 



CONTENTS 



5 



review of those results. 

This review aims at becoming a seed for a growing research field. It is organized as 
follows: in section[2]we will derive the stochastic Schrodinger and the master equation 
both in the Markovian and non-Markovian approximation. We will also discuss 
some of the issues in solving numerically the stochastic equation. This ingredient 
is fundamental if we want to discuss some applications to real systems. Section [3] will 
present some recent results on how to simplify the numerical investigation of many- 
body open quantum systems with techniques from the Density Functional Theory [13] . 
We will also discuss the possibility of closing the Kohn-Sham system as suggested 
recently [Tl] . 

Section |3] deals with some examples of application of the stochastic equation 
to real systems: we will discuss the case of spin thermal transport, Bose-Einstein 
condensation, thermal relaxation and the effect of electron energy dissipation on the 
ionic motion. 

2. Open Quantum systems 

The theory of open quantum systems has a long history that dates back to the 
beginning of quantum theory. The inclusion of the coupling between a quantum system 
and an external environment reached a certain degree of maturity with the pioneering 
works of Vernon and Feynman TB] together with Caldeira and Legget [TB] . The general 
idea is to derive an effective dynamics for the quantum system that takes into account 
the coupling with the environment without solving the equation of motion for the 
environment. The starting point has been the von Neumann equation for the density 
matrix, since the latter can be easily connected to the thermodynamical properties of 
the system pTl [TS] . Within the so-called master-equation formalism, an impressive 
number of results have been obtained, setting it to be the "de facto" standard for the 
theory of open systems. 

More or less in parallel as what happened in the classical theory with the Langevin 
and the Fokker-Plank equations, one can derive an effective dynamics for the "state" 
of the quantum subsystem, the so-called stochastic Schrodinger equation (SSE). The 
theory began with an attempt to, mimicking the situation of the Langevin theory of 
Brownian motion, introduce a stochastic term to describe the relaxation dynamics of 
an open quantum system [TTl [Tni HOI [5TJ [231 UHl HH IH] ■ It then found application into 
the theory of quantum optics where the environment is the electromagnetic radiation 
|26j . The stochastic equation was meant to reproduce the dynamics derived from the 
density matrix after some average was taken. Recently, Gaspard and Nagaoka [57] 
formalized the theory starting with the equation of motion for the environment, the 
system and their coupling. With some approximations, which are similar to those 
invoked for the derivation of the master equation, they arrived at a general expression 
for the stochastic Schrodinger equation [37] . The idea behind this approach is similar 
to the Gibbs' ensemble theory [TT] HH HH] : we build many replicas of the same system, 
each identified by a certain realization of the dynamics of the environment -assumed in 
thermal equilibrium- and let them evolve independently. To obtain physical quantities 
from this amount of information we perform averages over the many realizations of 
the micro-state of the environment. 



CONTENTS 6 
2.1. Stochastic Schrodinger Equation 

We begin by considering a subsystem described by the Hamiltonian Hs coupled to 
an external environment, given by Hb, through an interaction potential XW, the 
Schrodinger equation for the closed system reads (?i = 1 hereafter in this review) 

id\^{t)) = {Hs + HB + XW)\^{t))dt. (1) 

We depict a possible situation in figure [TJ where a system is coupled to three different 
baths. Each of the baths is characterized by its own thermodynamical micro-states as 
we will discuss in the following. 




Figure 1. The subsystem S is in contact with three external bath Bi, B2, and 
-B3 , that together form its "environment" . Each of the bath is characterized 
by its micro-states, and possibly by a global thermal equilibrium at given 
temperature. We allow for energy and momentum exchange within the system 
and the environment. 

This differential equation describes the exact dynamics of the closed system. 
However, the exact microscopic description of the dynamics of the macroscopic 
environment and its influence on the subsystem are in most cases not feasible. 
Consequently, this equation will serve as a starting point for the derivation of an 
equation of motion for the reduced wave function expressed in the Hilbert space of 
the subsystem. The following deduction is very much in the spirit of Gaspard and 
Nagaoka ^7\, who applied a so-called Feshbach projection-operator method ^51 [50] 
to the Schrodinger equation ([1} and derived a non-Markovian stochastic Schrodinger 
equation. 

By considering a complete and orthonormal basis for the environment, 

1b = ^\n){n\, HbIu) ^ en\n), {m\n) ^ 6„in, (2) 

n 

the total wave function can be expanded in this basis as \'^{t)) — J2n ® k^)- 

Due to the normalisation of the total wave function, 

1 = (*|*)=^(0"|0"), (3) 

n 

the coefficient wave functions are not normalised and the square of their norm 
can be interpreted as the probability that the environment is in the state |n). As a 
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consequence, the wave functions form a statistical ensemble which determines the 
total state of the combined system. In order to extract a typical representative of the 
ensemble, one defines the projection operators 

P=ls<E>\l){l\, Q=ls<E) J2\n){n\. (4) 

Hence, by applying the operator P to the total wave function, we extract the l-th 
coefficient wave function, PI'S) = \(pf) (E) Correspondingly, Ql^) contains the 
information of the other wave functions belonging to the ensemble. We want to point 
out that these projection operators satisfy 

p^^P^pIq^^q^q^^P + q^It, (5) 

and thus can be used for the Feshbach projection-operator method. Conveniently, this 
method is performed in the interaction picture 

id\^i{t)) = XW{t)\^'i{t))dt, (6) 

where the total wave function and the potential in this picture are given by 

\^i{t)) =e*^«*e*^^*|*(t)), 

W{t) = e*-«"«*e*-^s*W^e~'-^^*e~*-^^*. (7) 

The idea of the Feshbach projection method is to split the Schrodinger equation for 
the closed system into two equations. One contains the information about the time 
evolution of a typical representative of the ensemble, the other is a differential 

equation for Ql^*) and describes the time evolution of all the other coefficient wave 
functions. Thus, by solving the second equation and inserting its solution into the 
first, one obtains a closed differential equation for To follow this plan, we apply 

the projection operators (|4]) to the time-dependent Schrodinger equation (|6]) and this 
leads to 

idP\^i{t)) = XPw{t)P\'fi{t))dt 

+ XPw{t)Q\^j{t))dt, (8) 
idQ\^i{t)} = XQW{t)Q\'i!i{t))dt 

+ XQW{t)P\^i{t))dt. (9) 

The second expression is an inhomogeneous linear differential equation for 
and can be solved by the method of variation of constants 

Qi*/(<)) - umi^iio)) 

-ixf dt'U{t-t')QW{t')P\^i{t')), (10) 



where U{t) is the time-evolution operator of the corresponding homogeneous 
differential equation and thus obeys idU{t) = XQW {t)QU {t)dt. Inserting (flO)) into 
([5]) leads to a closed differential equation for 

idP|«'/(i)) = XPw{t)P\^'i{t))dt 

+ XPw{t)(^U{t)Q\'ii{0)) 

-iA^ dt'U{t-t')QWit')P\^i{t'))^dt. (11) 
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It is worth mentioning that until now no approximations have been made, so that (|lip 
describes the exact time evolution of the Z-th coefficient of the total wave function. 
Hence, (|TT]) can also be considered as a suitable starting point for a derivation 
of a stochastic Schrodinger equation beyond the weak coupling approximation. 
Unfortunately , this equation is as difficult to solve as the Schrodinger equation for the 
closed system ([T]) and thus we will consider a subsystem which is weakly coupled to 
the environment. To this end, we perform a perturbation expansion to second order 
in the coupling parameter A, 

idP\^i{t)) = \Pw{t)P\^!i{t))dt + \Pw{t)QYii{Qi))dt 
-i\^Pw{t) dt'{QW{t')Q\^i{Q)) 

+ QW{t')P\^i{t'))^dt + 0{X^), (12) 

where the time-evolution operator U (t) has been expanded to second order in 
A. Until now the derivation was quite generic, no restriction on the form of the 
interaction potential W or the Hamiltonians of the subsystem and environment has 
been imposed. In the following we assume the interaction potential to be of linear form, 
W = J2a ^aBa, in the operators Va and of the subsystem and the environment, 
respectively. These operators can always be redefined as Hermitian operators [27], 
thus we use = Va and = Ba in the following. If needed, this restriction can 
easily be lifted. 

By multiplying (|12p from the left with {l\ and assuming that {l\Ba{t)\l) vanishes 
it simplifies to 

id\cf>\{t)) = \f{t))dt 

tX^y^Vait) f dt'Vp{t') (13) 



a, 13 

4.1 



X {l\Ba{t)Bp{t')\lM{t'))dt. 
The forcing term 

1/(0) = A y Vait){{l\Bc.{t)\n) 

~ i\ j\t%{t'){l\Ba{t)Bp{t')\n)]\4>'}m (14) 

describes the influence of all the other bath modes on the /-th coefficient wave function 
and one sees that the initial conditions |(/)"(0)) enter here as an essential ingredient. 
By assuming that at i = the subsystem is in a pure state and the bath is in thermal 
equilibrium, the total density operator can be written as 

pt{o) = \m){m\^p''^ = 10(0)) (0(0)1 ® (15) 

where /?' is the inverse of the temperature and Zb = Tr^ e^^ . Nevertheless, we 
are interested in the initial wave function corresponding to this density operator. This 



§ This condition can be either fulfilled through a redefinition of the systems Hamiltonian or 
through the choice of the operators Ba ■ 
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can be established under the assumption that the initial condition for the total wave 
function is given by 



\^>m^\m)®Y. 



Zb 



(16) 



where {On} are independent random phases uniformly distributed over the interval 
[0, 27r]. Here, we want to point out the difference in the description of an open quantum 
system by wave-function or master equation methods: random phase factors have to 
be considered in the initial conditions. Due to this, the initial conditions can be 
written as 



|^"(o)) = (H*(o))-|0(o))i 



Zb 



. f' 



(17) 



where we have used the fact that all coefficient wave functions att — are proportional 
to the same state |0(O)) of the subsystem and thus they can be expressed in terms of 
the l-th. With the help of this, the forcing term ((Ti| can be simplified further 

i/(t))«A J2 Ut){i\Bc.it)[\^'m 



P 



A J2 V^it){l\Bc.{t)\nM{t) 



g--(e.-ei)g*(e 



(18) 



Here, we have used that the expression in the square brackets gives the time evolution 
of the l-th. coefficient wave function in the interaction picture up to second order in A 
for ([T5)) . Besides, we have included the stochastic noise in 



(19) 



which depends on a specific coefficient wave function. In order to eliminate this 
dependence, a thermal average is performed. 



e ~ 



rZs 



1 



iz^ 



(20) 



Additionally, one assumes that the expectation value of an operator from the bath 
for a typical eigenstate |^) is approximately equivalent to a thermal average of the 
temperature of the bath. 



{l\B^{t)Bfin) « Trs if^B^{t)Bt3{t') = Co,.p{t - t'). 



(21) 
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For a more complete analysis on the validity of this assumption see [31] |3H 1331 1311 
[35l l36] . In this expression we have defined the bath correlation function Cq^^, which 
describes the influence of the bath onto the system by tracing out the bath degrees of 
freedom in the dynamics. 

In addition, if the bath is large enough, 7q, {t) consists of a sum of many complex 
oscillating terms which leads to random Gaussian behaviour according to the central 
limit theorem. Hence, the noise is characterised by its mean value and its variance. 



7«(i) = 0, lo.{t)i0{t') - 0, 7^(i)7/3(t') = C^a,M* " (22) 



where the relations e^(^"+«'") = 0, e*(^""^"-) = 5nm and {l\Ba[t)\l) = have been 
used. We want to point out that the noise and the bath correlation function are not 
independent, more precisely, the covariance function of the noise is given by the bath 
correlation function. Collecting all the information, transforming back into the partial 
Schrodinger picture of the system and setting t ~ t ~ t' , (I13|) can be written as 

id\m) = Hs\<t>{t))dt + \Y,io.{t)v^m))dt 

a. 

-iX^Y.^-'l dTe-'"'^%C^Ar)m-T))dt. (23) 

Here we have suppressed the index I, since we assume this wave function is a "typical 
representative" of the dynamics of the system. This again corresponds to the Gibbs 
ensemble theory: with probability close to 1, we are sure that picking at random one 
of the coefficient wave function, it will evolve according to In this equation 

the change of the wave function at time t depends not only on the current state 
\4'(t)) but also on an integral over the whole history of the state in the interval [0,t]. 
This behaviour is called non-Markovian and thus (j23p is denominated non-Markovian 
stochastic Schrodinger equation (NMSSE). 

In (j23p the coupling of the bath to the subsystem is described in an approximate 
manner and enters in the NMSSE through the bath correlation function Ca,p{t) and 
the stochastic noises ^a{t) with the properties ((22)) . Thus, all the information about 
the time evolution of the bath and its coupling to the subsystem is included in the 
bath correlation function. In most quantum optics cases the dependence on the 
past of the wave function can be neglected. This is due to the fact that the bath 
correlation function decays rapidly to zero on a time scale on which the system wave 
function does not vary significantly. For convenience, one neglects the non-Markovian 
behaviour by approximating the time dependence of the bath correlation function by 
a ^-function, Ca,p{t — t') « Da^p5{t — t'), known under the name of a 5-correlated 
bath approximation. As a result, the NMSSE reduces to 



tdim = 



X \(j){t))dt. (24) 



where 7a (t) are white-noise processes with 7a (i) = and 7a(07^i(^') = Dci,i3S{t — 
t'). With the help of a unitary transformation U^^s, that diagonalizes -Da,/3 with 
corresponding eigenvalues da, (|24ll can be written in an Ito differential form 



dm)) = [[ - iHs - \Y.^^o.Sa\dt + Y.^o.dW^^\m), (25) 



CONTENTS 



11 



where the new set of bath operators is given by Sa = ^V^a^p Ua.fjVf}. Furthermore, 

the stochastic processes are included in dWa = ^"i/Vd^J^j ^] alj'^^ ^^'^ ^''^'^ '-^'^ show 
that these satisfy the properties 

dW2 = 0, dWadW*p = Sa^pdt. (26) 

We will call Markovian stochastic Schrodinger equation (MSSE) and we want 
to point out that this equation does not follow standard rules of calculus. The state 
\4>{t)) is a stochastic function and its time derivation is not defined at any instant of 
time. In addition, the differential noise dW scales on average as Vdi and thus can be 
interpreted as a differential increment of an underlying Wiener process. As a result, 
this differential equation is not tractable with standard calculus and the rules have to 
be modifies according to the ltd calculus. Since this will bring us too far from the scope 
of this review, here we will only state the important results of the Ito calculus needed 
in the following sections. For a more complete treatment of the Ito formalism one 
can consult the vast literature on the subject, here we just point out a few standard 
references [371 133 ■ An important result from this calculus is the Ito chain rule [3] [35] , 

d{\4>){i,\) = (7^1 + + (27) 

where and if] are two states evolving according to the MSSE In addition, the 

rules for the Ito differentials 

dWadWj = Sa,pdt, dWcdWfs = 0, dW^dt = 0, (28) 

should be kept in mind. In the following section we will apply these Ito rules in order 
to derive the master equation that corresponds to the MSSE. In the case of the non- 
Markovian SSE, where one is confronted with non-white noises, Ito calculus cannot 
be applied and one has to find another approach to derive the corresponding master 
equation. 



2.2. Density Matrix formalism 

In the previous section we have discussed the dynamics of a quantum mechanical 
system coupled to an external environment from the point of view of what we have 
identified as the "state" of the system, However, for historical and practical 

reasons this is not the standard starting point. It has been easier, as we will discuss 
in a moment, to start from the dynamics of the reduced density matrix or statistical 
operator of the system, defined as 



p=\4'){<P\. (29) 

To obtain the equation of motion for the density matrix that corresponds to the 
Markovian SSE (I^S)) . we can start from pS)) and calculate the differential 



dp = = {d\q^)){^\ + + (30) 

Unlike in normal calculus, one also has to keep the term d\4>)d{(j)\ as it contributes 
on average to first order in dt. In the Markovian case, where one has to deal with 
white- noise processes, we can apply the Ito rules (l?fl) and (1^5)) which lead us to 



dp -■ 



iHsdt + ^[ - \siSadt + SadWa^ 



(31) 
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+ P 



J2 So^dWapSldW* + 0{df 



i[Hs,p]dt - i ^ IsiSaP + pSiS^ - 2S^pSi\dt + 0{dr 



This well-known Lindblad master equation is the most general type of a Markovian 
master equation which is known to preserve not only the norm but also positivity and 
hermiticity [JO] . This means that the Markovian SSE describes on average an open 
system dynamics which coincide with the Lindblad dynamics. Here, we have assumed 
that the Hamiltonian of the system is non-stochastic. However, if the Hamiltonian is 
stochastic, one has to deal with an ensemble of Hamiltonians and the corresponding 
equation of motion will most likely differ from the Lindblad master equation. We will 
discuss this issue in depth in section 12.51 where a gas of interacting bosons will be 
considered. 

As mentioned before, Ito stochastic calculus is only applicable for white-noise 
processes. However, in the non-Markovian SSE (1^^ one encounters with coloured 
noise which is not (5-correlated. A stochastic calculus for coloured noise is not so 
extensively investigated as the Ito calculus. Hence, one has to find an alternative 
procedure of deriving an equation of motion for the density operator that corresponds 
to the NMSSE. To this end, a comparable master equation for the NMSSE up to 
fourth order in A can be derived by performing a perturbative expansion in A to arrive 
at HZ] 

= -^[Hsrp{t)\ (32) 



dt 



A^ /* dr ^ Clp{T)V^p{t)e''"''Vpe^^'^ 
A^ /* dT V Ca.Ar)e-'"''%e'"'^ p{t)Vc. 
A^ /' dT Cc.Ar)Vae-'"'^V(ie'"'^p{t) 



In the following, we will call this equation non-Markovian master equation as it 
corresponds to the NMSSE. Also here, we have assumed that the Hamiltonian is 
non-stochastic. In the limit i — oo in the history term, the non-Markovian master 
equation simplifies to 

^=-^[Hs,m] (33) 
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where 

/>oo 

Ka^X'y2 dTCa^,f,{T)e-'"^^V^e'"'\ (34) 

This time-local evolution is the well-known Redfield master equation [JT] and has 
been widely applied to systems where the dynamics of the environment is much faster 
than the system dynamics. This equation in general provides information on the long 
time relaxation dynamics of the density matrix. In conclusion, we have seen that the 
Markovian SSE describes on average a dynamics that corresponds to the one obtained 
from a Lindblad master equation and the non-Markovian SSE corresponds to the 
master equation if the Hamiltonian is not stochastic. Hence, the SSEs can also 
be seen as an unraveling of the corresponding master equation, i.e., a quick and dirty 
way to obtain the solution of a master equation, especially when the systems size 
grows one expect a better scaling behaviour of the SSEs. 



2.3. Equation of motion for observables 

Given the equation of motion for the density matrix and the SSE, either in the non- 
Markovian or in the Markovian approximation, the next step is to derive a general 
equation of motion for the expectation value of any observable. By definition, given a 
physical observable P, described by the operator P, its expectation value is given by 



ip^^mEm^TrsmP) (35) 

(mm 



N 

I 

lim 



1 {Mt)\P\Mt)) 



where, again, the (. . .) represents the standard quantum mechanical average and 
the average over the noise as we have introduced them in section [2. II To understand 
this last process, we can think of pSI) as having build a certain number of replicas of 
the system, labelled by the index i, (j)i{t), see figure [H Each replica evolves according 
to the SSE with a given realization of the noise. Averaging over the noise therefore 
corresponds to summing up all the weighted contributions to the given observable 
coming from each replica. 




Figure 2. In evolving tiie SSE we build many replicas of the same system and let 
them follow their own time evolution. In building up the observables, we average 
over the contribution to the observable of each replica, as in II35II . 

It should be evident from psp why the density matrix got a relevant amount of 
attention: if one were able to solve the equation of motion for the density with the 
proper initial conditions, the expectation value of any physical observable would come 
essentially at no extra cost. From it is clear that the process of averaging over 
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the many realizations of the stochastic behavior that is driving the dynamics of \^) 
is a crucial part of the definition of the density matrix. It is this average that builds 
up the correlation between the different trajectories of the stochastic dynamics. The 
same average step makes the SSE difficult to handle: even if one were able to write 
down an explicit solution to the SSE, the physical quantities will appear only after the 
average over many realizations of the noise is taken. For the density matrix, this step 
is embedded in the definition, and is thus carried out before the dynamics is solved, 
obtaining therefore an equation for the averaged quantities. 

We want to derive an equation of motion for (P). To this end, we consider the case 
of the Lindblad equation or the MSSE (P5|) . The generalization to the non-Markovian 
dynamics is lengthly involving the expansion of the equation of motion in terms of the 
coupling parameter A [27] . 

Considering the case of a single Wiener process, dW , we have by using ltd calculus 

d{P) = + (0|P(d|0)) + 

-i\P, H]-]-( S'^SP + PS^S] )dt+( S^PS) dWdW 



2 

+ {S^P + PS)dW + 0{dWdt). (36) 

The equation of motion for the ensemble-averaged expectation value is obtained 
immediately from ([55)) . 



dt{P)^ -^{[P,H]) 

- i ((S'tS'P) + (PS'tS') - 2(S'tPS')) (37) 



P,H 

- i (^{S^SP) + (PSIS) - 2{S^PS)^ . (38) 

One could arrive at the same equation by using the master equation for the density 
matrix. It is immediately seen that the coupling to an external environment, modifies 
the dynamics of the observable if P and S do not commute, i.e., [P,S] ^ 0. The 
second term in the right hand side of (155)) describes the dynamics of the expectation 
value in the presence of the relaxation induced by the environment. A steady state for 
any observable will be reach when a competition between the two terms is established 
so that to make the right hand side of (|38l) vanishing. 

Of importance for the investigation of the dynamics of open quantum systems, is 
the time evolution of the single-particle density, i.e., the expectation value of the single- 
particle density operator. In standard quantum mechanics, i.e., for closed systems, 
the equation of motion for the single-particle density is the continuity equation which 
essentially embodies the important principle of particle conservation. Namely, if n(r, t) 
is the single-particle density and j{r,t) the single-particle current density, we have 

dtn{r,t) = -V ■j{r,t). (39) 

When the system is open, the equation of motion for the single-particle density appears 
more complex. We have from (|55)) that 

dMr^ = -V • J (r, t) - i ((S-t^^) + (hS^'S) - 2{S''hS)^ . (40) 

It may appear that this equation breaks the important physical principle of particle 
conservation, since the second term in (HHj) can hardly be recognized as a divergence 
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of any current. Improbable as it might seem, Gebauer and Car |42| were able to show 
that, at least in the Markov approximation, the right hand side of (j40p can be written 
as the divergence of a new current density, j + jc , where jc describes the transfer of 
momentum between the system and the environment. 

It will be important when we will discuss in the following section the development 
of a density functional approach to open quantum systems, so let us discuss something 
that might appear as a trivial point. From a mathematical point of view, the 
continuity equation, given the initial condition and the current density j(r, t), uniquely 
determines the particle density everywhere and at any time. This is because a simple 
integration of (15^ gives 

n{r,t) ^n{r,t = 0)- [ dt'V ■ j{r,t'). (41) 
Jo 

The opposite is not true. Given n{r,t) and the initial conditions, the current 
density j(r, t) is not uniquely determined. Indeed, given j{r, t) and j{r, + V x v{r, t), 
where u(r, t) is an arbitrary vector in space and time, these two currents solve the 
continuity equation with the same density n{r,t). For this reason, one usually says 
that the particle density only determines the longitudinal current-density, while the 
transversal part is left unknown. Indeed, only the longitudinal current is responsible 
for the transfer of particle across any given surface. 

The situation appears more complex if one starts with the continuity equation 
for open systems. Even if one would consider the total current j + jc, the density is 
not uniquely determined. This is because the current jc is determined by the density 
and the current density themselves, therefore making (j40[) a non-linear equation whose 
solution appears non-trivial. The situation has not been clarified and in the following 
we will assume that, similarly to what happens with the continuity equation, (|40p 
uniquely determines the single-particle density. 

With logical steps similar to those who led us at the equation of motion for the 
single-particle density we can derive the equation of motion for the ensemble- 
averaged current density for a system of interacting identical particles of mass m, in 
the presence of an external vector potential Af,xt, 



dtj{r,t) = -^^dtA,,t{r,t) - '-^^ x V x A,,t{r,t) 



{g{r,t)) (42) 

lit 

where we have defined 

Gir, t) = S^j{r, t)S - ij(r, t)S^S ~ i^t^i(r, t\ (43) 

-^('^' = - ^ '5(r - f,)V, t/„,t (f, - fj) + mV • a(r, t) (44) 

with the stress tensor a(r^ t) given by 

Hr,t) ^ -^Y.i^^^i^i^^^'^ - ^k)}}. (45) 

In these equations Vj contains the derivatives with respect to the coordinates of the 
j-th particle, i.e., in 3D Vj — {dxj,dy-,dzj), and the current operator is defined as 

= jEi'^l^-^.),^^^} (46) 
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with 

Pi + eAext{fi,t) 

Vi = , (47) 

m 

the velocity operator of particle i, and the symbol {A, B} = {AB + BA) is the anti- 
commutator of any two operators A and B. 

The first two terms on the right hand side of (UU describe the effect of the applied 
electromagnetic field on the dynamics of the many-particle system; the third is due to 
particle-particle interactions while the last one is the "force" density exerted by the 
bath on the system. This last term is responsible for the momentum transfer between 
the quantum- mechanical system and the environment. It should be pointed out that 
the first three terms on the right hand side of (|^^ are present also in the standard 
equation of motion for the single-particle current density [33] . Equation (H^ can be 
seen as the equation for the "forces" acting on a single particle. 



2.4- Solving the Stochastic Schrddinger Equation 

For a single realization of the noise, the numerical solution of the SSE is similar to 
the integration of the Schrodinger equation [57, , and one can use some of the 
standard techniques for the numerical integration of partial differential equations [46] . 
However, the noise effectively reduces the stability and efficiency of the numerical 
algorithms |44) . For this reason, higher order techniques, like the standard fourth 
order Runge-Kutta, do not offer the usual improvement as in the integration of the 
standard dynamical equations. Indeed, the usual parameter used to identify this 
efficiency is the scaling of the numerical error, e, in terms of the time step of the 
numerical integration. At. Typically, we have a power law scaling, i.e., 

e cx At''. (48) 

For example, for ordinary differential equations, the second order Runge-Kutta has 
V = 2, while the fourth order Runge-Kutta has v = A. For Markovian systems, it is 
found quite surprisingly that for any standard technique, e.g. the second and fourth 
order Runge-Kutta methods, we have the same v [4 7) . Indeed, it has been shown 
that, in general, any method which involves only Wiener processes, will have the same 
scaling [17] as the simpler Euler or Heun scheme (second order Runge-Kutta). High 
order schemes which improve the scaling with the integration time step can be found, 
and have been proposed -see for example [HJ US] and references therein- which involve 
the evaluation of a class of stochastic processes, beyond the simple Wiener processes, 
in building the numerical approximation to the exact solution of the SSE. These terms 
appear in the Taylor-like expansion of the solution to high orders [331 |3H] • On the other 
hand, to the best of our knowledge, for non-Markovian SSE there are not clear results 
and a complete analysis of the error scaling is missing. The reader should also be 
warned of some attempts at building high-order schemes to the solution of Markovian 
SSE, like in [49] starting from the standard Runge-Kutta schemes. These are based 
on the assumption AVKAVF — At, which is however valid only when the average over 
many realizations of the noise is performed. Although it seems these schemes improves 
the convergence, their validity is doubtful [SOI and most likely restricted to a small 
class of equations. 

In this review, we will not discuss in details any algorithm of numerical solution of 
the SSE. A nice step-by-step tutorial on the solution of stochastic differential equations 
can be found in [3S], while more advanced techniques can be found in [331 [SO]- Here, 
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we would like just to add that we can find a non-Unear SSE which preserves, at each 
time step, the norm of the wave-function and which reproduces the same physical 
quantities as the linear SSE [51]. It is our experience that using this non-linear SSE 
usually improves the stability of the integration algorithm. 

When dealing with the non-Markovian SSE, we face the problem of simulating 
the correlated colored noise ^{t) in ^F^. Many solutions to this problem have been 
proposed: they do differ on the amount of information we need at our disposal about 
the function Ca^pit). For example, for simple correlation functions the algorithm 
proposed in [52] is probably the most efficient, but it does require the knowledge of 
a first-order differential equation, whose Ca,p(t) is a solution. This works well for 
Ca,i3{t) oc e~l*l [S2]- For more complex cases, this approach can not be used, and 
we can revert to the solution proposed by Rice [53], and recently revisited [53]. This 
algorithm might not be the most efficient [551 ISS] ; but is based on the sole knowledge 
of the Fourier Transform of Ca^p{t), the power spectrum, practically the minimal 
amount of information we need. Moreover, this Fourier Transform can be known 
analytically for some models: this is due to the fact that we might have access to, 
or we can approximate it to a certain degree, the power spectrum of the bath. A 
variant of this algorithm has been recently proposed in |57) for the case we know the 
function Ca,p{t) at each instance of time, and we do not want to store the full Fourier 
Transform for computational purposes. However, this algorithm does have a large 
numerical overhead that makes its application to the SSE too expensive [58] . 



2.5. Difference between the Stochastic Schrddinger equation and the master equation 

We want to discuss the equivalence between the master equation for the reduced 
density matrix and the SSE. In particular, we want to show one case of general interest 
in which the two formalisms are giving different results. We argue that this originates 
from the different ways they deal with interaction. 

Let us consider the Hamiltonian of a one dimensional boson system (in second 
quantization) which, when the bath is not present, reads 

H = H0+ i?™* ^ Jdx ^Hx) (-^^ + V^ext(x)) ^(X) 

dxdx' i;'' {x)'tl;Hx')Uint{x - a;')V(a;')V'(a;),(49) 

where ipix) destroys a boson at position x, Vext{x) is a confining potential, and 
Uint{x — x') is the boson-boson interaction potential. We consider here the case of a 
diluted gas of Alkali atoms, since this is a system of great interest [591 IMl EI). For 
these gases the interaction potential can be substituted with the contact potential, 
i.e., 

U^nt{x ~ x') = go{N - l)5{x - x') = ~g5{x - x') (50) 

where go = Anh^a/m is determined by the scattering length a of the boson-boson 
collision in the dilute gas and N is the total number of bosons in the trap, so that 

M\=im- 

In what corresponds to a Hartree approximation for the boson wave function, we 
can go from the equation of motion for the annihilation operators to the equation of 
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motion for single-particle wave function ^'(x, t), when the external bath is not coupled 
to the boson gas, 



'2 



+.gn(a;,i)*(x,t) (51) 



where n{x,t) = |vl'(a;,t)p is the single-particle density of the boson gas [62j [63] . 
Equation ([5T|) (and its generalization to two and three dimensions) correctly describes 
the physical properties and especially the dynamics of a Bose-Einstein condensate 
[Sni [SOI [HI] ■ That is, ^'(a:, t) describes the dynamics of the ground state of the system, 
when the temperature of the gas is well below the critical temperature for the Bose- 
Einstein condensation [18l [171 EI] ■ 

A harmonic confinement is routinely generated in the magneto-optical traps used 
in the realization of the condensation in Alkali gases [SHJ [501 III] i so we choose 

Vext{x) = ]^mu)lx^. (52) 

When the boson system is coupled to the external environment, we assume that 
the Hamiltonian is not affected by this coupling, i.e., we neglect the Lamb shift, and 
the system evolves according to the Markovian SSE 

1 <e 1 

h -1 

2m dx"^ 2 



d^{x, t) = -i[ -7^73^ + ^mulx^ + gn{x, t) ) ^{x, t)dt 
1 



^S''S-^{x,t)dt + S-^(x,t)dW, (53) 

where we have assumed that the coupling with the environment is given by a single 
operator S. Here, dW is a standard Wiener process with properties (pS)) . 

For convenience, we rescale this equation in terms of the physical quantities wq, 
Xf) = 1 / ^mujQ , and g = g/xg to arrive at 

Xo__f_ ^ x^ ^ g 
2 dx^ 2xq LdQ 



d^>{x,t) = -iujol — ^-r^ + TT^ + —n{x,t)xo '^{x,t)dt 
\ 2 dx^ Ix 

1 



^S'^S^{x,t)dt + S^{x,t)dW. (54) 

We begin with considering the case of non-interacting bosons, i.e., we set 3 = 0. In 
this case, the Hamiltonian admits a natural complete basis, the set of Hermite-Gauss 
wave- functions 

^,(x) = , M x/x,)e--''^<, (55) 

where Jij^x/xo) are the standard Hermite polynomials |64| . If we expand the wave- 
function \E'(a;,t) = J2j 0'ji't)'Pjix), and make use of the orthonormality properties of 
the Hermite-Gauss wave functions, we obtain the (stochastic) dynamical equation for 
the coefficients aj, 

da, = J2 (H°,a, + ^{S^S),,a,) dt + dWj2 S^,a, (56) 
where iJ"- = (j + l/2)a;o(5y . 
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Together with (|56p we can study the dynamics of the density matrix via the 
quantum master equation pip , which in the same representation as (|56p reads 



k 

^ (^S^kPkk'Slj - -^Sji^Skk'pk'j - ■^PikSl^.,Sk'j^ ■ (57) 



k 



k,k' 

The connection between (I56p and (|57p is estabhshed by the identity pij = a*aj which 
is vahd for any pair of indexes i and j . It was shown |65j that when the interaction is 
set to vanish, the two formahsms yield the same result in the hmit of a large number 
of realizations. 

When one turns on the particle- particle interaction Uint, this corresponds 
to adding to the free Hamiltonian an interaction part, iJ™* that in the basis 
representation of the Gauss-Hermite polynomials reads 

^5* = E^^^;m4«, (58) 

k.q 

where Fi,j.k,q is the fourth-rank tensor defined as 

9 



2-(l+j+k+q)/2 



00 

X / dx Hi{x)Hj{x)Hk{x)Hq{x)e-'^''\ (59) 



A long but straightforward calculation, worked out in full detail in |66) . brings us to 
an expression for Fi^j-^k.q in terms of Euler gamma functions and a hypergeometric 
function [661 164) . It can be shown that the hypergeometric function reduces to the 
summation of a few - at most min(i, j) - terms [65) . In the case of the density matrix 
approach the interaction Hamiltonian is immediately written as 

k,q 

In solving the dynamics of the system described either by the SSE ([51]) or the 
master equation ([57| when the interaction is present, we have assumed that at any 
instance of time the bath operator brings the system towards the instantaneous ground 
state of the total interacting Hamiltonian H^^ + H"j* . Formally, in the basis set that 
diagonalizes the total Hamiltonian, this corresponds to choosing 



fO 1 1 ... 1\ 
... 



(61) 



V ... / 

for the bath operator. In addition, the interaction potential (and hence the total 
Hamiltonian), being defined in terms of the instantaneous density, is stochastic, 
namely it is different for each element of the ensemble. While we take this into 
account explicitly in the SSE ([M)) . in the master equation ([57| we must consider the 
interaction Hamiltonian averaged over all realizations. 

In figure 13] we plot the occupation probability Pj{t) of the state j for the first 3 
energy levels of the free Hamiltonian {pj{t) = \aj{t)\'^ from the SSE or Pj{t) = Pjj{t) 
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Figure 3. Probability of occupation of the first 3 energy levels of the interacting 
Hamiltonian versus time calculated via the SSE (solid lines) averaged over 100 
independent runs, and via the equation of motion for the density matrix (dashed 
lines). In the inset we compare the equilibrium density (black, dashed line) with 
the one obtained from the Thomas-Fermi approximation to the ground state 
(orange, solid line). The initial condition is chosen by putting the boson with 
equal probability in all the considered states. 



from the density matrix). We have assumed an interaction strength g/ujQ — 5, 
considered the first 20 energy levels, used a time step woAi — 60/2^^ and we have 
performed 100 independent runs of the SSE. The initial condition is such that all the 
energy levels are occupied with the same probability. From the figure it is evident that 
the system reaches, in the long time, the same steady state, but it is also clear that 
the state calculated with the SSE relaxes slower to this steady state than the state 
obtained from the density matrix equation. This is a spurious effect in the density 
matrix dynamics where the average density defines the interaction potential without 
account for the fiuctuations of the state, and hence of the stochastic Hamiltonian. 

We have also tested that the steady state reached during the dynamics is 
consistent with the eigenstates of the Gross-Pitaevskii equation [51] [ST]. In particular, 
the ground state of the interacting system, when the interaction is strong, can be 
obtained by neglecting the kinetic contribution to the Hamiltonian. In this case, a 
good approximation to the ground state density reads |61| 

^-0 X = °—e U - \ Imujlx^) 

+ terms proportional to l/g^, (62) 

where /i is the chemical potential, i.e., in this case the energy of the ground state, and 
e(x) = if a; < and e{x) = 1 if a; > 0. 

In the inset of figure [3| we plot the density obtained at tujQ — 60 from the SSE 
(black, dashed line) and the density obtained from the approximation (1621) (orange, 
solid line). Notice that the value of the parameters g and ^ have been obtained from 
the best fit with the numerics: indeed one can show that the approximation (|62p 
is exact in the limit of very large interaction, |67[ 161] which is not reached in our 
calculations. 

We have therefore shown that the master equation and the SSE, for a system 
where the Hamiltonian is stochastic, i.e., it does contain terms that depend on the 
wave function, could lead to different results. This is expected since in the master 
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equation we replace the effective stocliastic term with an "average" contribution, 
without having any control on it. For example, in the case discussed the approximation 
corresponds to replacing terms like n(x,t)n{x' ,t) with n{x,t) x n{x',t). In the SSE, 
this approximation is not performed, since the average over the stochastic noise is 
done after the full evolution is calculated. 

3. Stochastic Density Functional theory 

As we will see for example in section 14.11 one severe limitation of the applicability 
of the theory is the presence of interaction between the particles of the system. This 
is nothing new since similar problems arise also in the case of closed systems where 
the presence of interaction makes obtaining the exact solution of the dynamics of the 
system an almost impractical task. However, in the theory of many-particle systems, 
whenever we are only interested in certain quantities, it has been found that the 
complexity of the problem can be greatly reduced, up to the level of making it treatable 
with numerical techniques. These results go under the name of Density Functional 
Theory (DFT) [551 EH [13] ■ Nowadays, there are different flavors of what one can call 
Density Functional Theories, and the interested reader can find better introductions 
to that field for example in |13) . For our purposes, it is sufficient to say that, if we 
focus on the dynamics of one observable, let us say the single particle density, then 
we can obtain this dynamics by investigating a system of non-interacting particles, 
called Kohn-Sham (KS) system in the presence of a certain external potential [701171) . 
This potential is built to give exactly the dynamics of the single-particle density of the 
real system fl3l . However, it is important to stress that any other quantity calculated 
within the KS system and which cannot be expressed in terms of the single-particle 
density alone, cannot, in general, be trusted. 

For the sake of this review, of even more interest is the case in which the current 
density is the quantity we would like to obtain. For this case, it has been shown that 
one can build a system of non-interacting quasi particles that, in the presence of an 
appropriate vector potential, is able to mimic the dynamics of the current density of 
the real system |72[ I73| . Again, one must add the usual caveat that in principle only 
the current density and all the physical quantity that can extracted from it with simple 
operations (like for example the single-particle density) have a physical meaning. All 
the others bear little resemblance to the same quantity of the real system, although, 
in some cases one can realize that some physical information can still be extracted 
from the dynamics of the Kohn-Sham system 

The advantage of this formulation is manyfold and so far reaching that W. Kohn 
was awarded the Nobel Prize in Chemistry for the initial formulation of the theory 
which was dealing with the properties of the ground state. Indeed, it essentially paves 
the way to the numerical investigation of complex many-body systems. It would be of 
great interest if one were able to bring together the formulation of the open quantum 
systems we have discussed so far and the KS theory of many-body systems. This 
result has been achieved recently by various groups [751 HSl IHSl [13] with different 
formulations of the theory. Also in this case, the parallel between the density matrix 
formalism and the SSE appears evident and not surprisingly in certain cases one can 
reproduce the result of one formalism into the other. However, as we will point out 
more clearly in the following, a KS formulation in terms of the density matrix of open 
quantum systems appears flawed |65| since at the very basic level the KS Hamiltonian 
does contain non-linear terms of the state of the system thus making the derivation of 
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a closed equation of motion for the density matrix problematic as we have discussed 
in section [2] 

3.1. Stochastic Time- Dependent (Current-) Density Functional Theory 

Let us begin by revising, without proof, the theorems of time-dependent density 
functional theory. We consider N non-relativistic particles of charge e, mutually 
interacting via the Coulomb interaction, in a time-dependent external potential, either 
a vector or scalar potential. We consider the general Hamiltonian 

. ^\p^ + eA,,t{h,t)Y N 

^ = E^ ^ — - + yeAh.t) + Y.u{h-i-j) (63) 

where Aexti^i t) is the external vector potential, VextiT, t) the external scalar potential, 
and U (r) is the particle-particle interaction potential. The many-body system evolves 
according to the time-dependent Schrodinger equation, 

i'^^Himi (64) 

When there is only a scalar potential Vext, we have the following theorems of 
Time-Dependent Density Functional Theory fTDDFT) [7P1 [7T1 17 ^ [771 : 

Theorem 1 (Runge-Gross [70]). We consider N non-relativistic electrons, mutually 
interacting via the Coulomb repulsion, in a time- dependent external potential. Two 
single-particle densities n(r,t) and n'(r,t) evolving from a common initial state 
\^{t = 0)) under the influence of two potentials Vext{i',t) and Vg^f{r,t) (both Taylor 
expandable about the initial time 0) eventually differ if the potentials differ by more 
than a purely time- dependent (r-independent) function: Vextir,t) — V^^^{r,t) ^ c{t). 
Under these conditions, there is a one-to-one mapping between densities and potentials, 
which implies that the potential is a functional of the density. 

The theorem is similar to the one that Hohenberg and Kohn |B8| proved in 1964 
and that allowed, a year later, Kohn and Sham [69^ to formulate the DPT. The key 
point being that one can obtain physical information about a many-body system by 
investigating the behaviour of a simpler system where particle-particle interaction is 
turned off. 

The initial formulation of the Runge-Gross theorem was not satisfactory since it 
was based on the existence of some action that in the dynamical case did not respect 
causality [71]. Some years later, van Leeuwen was able to extend the theorem and 
prove that there is no need to define an action from where the equation of motion 
can be derived. Recently, Vignale has shown that the initial formulation can be 
made consistent by considering carefully the boundary and initial conditions for the 
dynamics [75] . 

Theorem 2 (van Leeuwen (IV). Let H{t) and H'{t) be two Hamiltonians of the form 
I163\) containing not only two different time- dependent local potentials VextifTt) and 
V^^^{r,t) but also two different particle-particle interactions U and U' . Let n{r,t) be 
the density that evolves from the initial state \^{t = 0)) under the Hamiltonian H{t), 
and let {"^'(t = 0)) be another initial state with the same density and the same value of 
the density gradient. Then the time- dependent density n{r,t) uniquely determines, up 
to a time- dependent constant, the potential V^^f{r,t) that yields n{r,t) starting from 
\^'(t = 0)) and evolving under H' . 



CONTENTS 



23 



It can be easily seen that the van Leeuwen's theorem include the Runge-Gross 
results as a simple corollary, by considering two systems evolving with the same 
interaction and with two external local potentials. 

The case for which a vector potential is present has been investigated first by 
Ghosh and Dhara [75] and later Vignale proposed a theorem along the proof of van 
Leeuwen. Here for simplicity we report only the second theorem, again the result by 
Ghosh and Dhara follows as a simple corollary. Vignale's theorem for Time-Dependent 
Current-Density Functional Theory (TDCDFT) states: 

Theorem 3 (Vignale [72 ). Consider a many-particle system described by the time- 
dependent Hamiltonian 

H{t)=Y^ ■^(p^+A,,t{h.t)f+Ve^t{f^,t) + ^ ?7(f, - ) , (65) 



where Vext{fit) and Aext{r,t) are given external scalar and vector potentials, which are 
analytic functions of time in a neighborhood oft^O, and U{ri, rj) is a translationally 
invariant two-particle interaction. Let n{r,t) and j{r,t) be the particle density and 
current density that evolve under H from a given initial state |5'(0)). Then, under 
reasonable assumptions, the same density and current density can be obtained from 
another many-particle system, with Hamiltonian 

^'W=E ^^{p^+A'if,,t)y + V'{h,t) +Y,U'[h-fj), (66) 

i L ^ i<j 

starting from an initial state |^'(0)) which yields the same density and current density 
as |5'(0)) at time t — 0. The potentials V'{r,t) and A'{r,t) are uniquely determined 
by Vei^t(r, t) and A^xtifit), 1^(0)), and |5''(0)) , up to gauge transformations of the 
form 

V'ir,t)^nr,t)-'-^, 

A'{r,t) ^ A'{r,t) +VA{r,t), (67) 

where A is an arbitrary regular function of r andt, which satisfies the initial condition 
A{t = 0) = 0. 

TDDFT and TDCDFT have been quite successful in describing the dynamics of 
closed quantum systems, especially whenever one needed to calculate the response of 
the system to an external perturbation [13]. These theories allow to go beyond the 
actual state-of-the-art in many fields, e.g., in electron transport in nanoscale devices 
where it does predict novel and interesting results [73] . The secret of this success is in 
the fact that the theory allows for treating complex system with relative ease, since 
only the dynamics of a closed non-interacting many-body system is computed. 

It appears therefore natural to try to generalize the theorem of TDDFT and 
TDCDFT to the case of open quantum systems. This generalization has been achieved 
by using both the reduced density matrix and the SSE formalisms [75] [751 |SS] . Let 
us assume that the quantum-mechanical system described by the Hamiltonian (1631) 
is coupled, via given many-body operators, to one or many external baths that can 
exchange energy and momentum with the system. If we assume that the dynamics of 
each bath is described by a series of independent memory-less processes, the dynamics 
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of the system is governed by the stochastic Schrodinger equation in the Markov 
approximation, 



+ Y,dW^it)Sc.\^{t)) 



where {Sa\ describe the couphng of the system with the environment. In 1681 we 
could have time-dependent Hamiltonian and bath operators. This would not change 
the following discussion. 

As discussed in section 12.31 we here assume that the knowledge of the current 
density is sufficient to obtain the single-particle density, and moreover this solution 
is unique, i.e., it depends solely on the single-particle current density and the initial 
conditions. If this is not the case, the proof of the following theorem, as we report it, 
does not hold and we have to revert to a more cumbersome yet similar proof where 
we need to assume that the ensemble-averaged single-particle density is difFerentiable 
at all order in time. 

We have the following result [75] : 

Theorem 4. Consider a many-particle system described by the dynamics in I168\) with 

the many-body Hamiltonian given by I163\) . Let n{r,t) and j{r,t) be the ensemble- 
averaged single-particle density and current density, respectively, with dynamics 
determined by the external vector potential Af.xt{r,t) and bath operators {So}- 
Under reasonable physical assumptions, given an initial condition |^(0)) and the 
hath operators {Sa\, another potential A'{r,t) which gives the same ensemble- 
averaged current density must necessarily coincide, up to a gauge transformation, with 
Aext{r,t). 

A theorem of similar content for the one-to-one correspondence between the 
single-particle density and the scalar potential has been discussed in [75] . The starting 
point of those Authors has been the equation of motion for the density matrix of the 
system, described in the Markov approximation. In [TU [50] the theorem by Burke et 
al. [75] and the previous one have been extended to the non-Markovian dynamics. 
The proofs of these theorems follow essentially the same logical steps as the ones we 
will present in the following. 

Proof: We follow a line of reasoning common to the proofs of the van Leeuwen 
and Vignale theorems [TT] [72]. Recently, a more elegant proof of Runge- Gross and 
van Leeuwen's theorems has been proposed [77], but the application of this method 
to the equation of motion for the current-density appears not clear. Let us continue 
by assuming that the current density j (r, t) is obtained from another many-particle 
system with Hamiltonian 

H'W - E ^ + o E - ' (69) 



2m 2 



evolving from an initial state |5''(0)) and following the stochastic Schrodinger 
equation (|68l) with the same bath operators Sa- By assumption, |vl''(0)) gives -in 
the primed system- the same initial current and particle densities as in the unprimed 
system. The equation of motion for the ensemble-averaged current density is (|42p . 
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Equations similar to (j45p - (|T7l) are written for the system with the vector potential 
A'{r^t). Similar force terms F' and Q' appear in these equations. F' and Q' differ 
from the forces in the unprimed system, since the initial state, the external vector 
potentials, and the velocity v are different. By assumption, the ensemble-averaged 
current and particle densities are the same, therefore 



V X A'{r,t) 



{r{r,t)) 



Taking the difference of and (1701) we arrive at 



n(r, t)dtAA{r, t) = j(r, t) x V x AA(r, t) + AQ(r, t) 



(70) 



(71) 



where AA(r, i) ee A'(r,t) - Aext{r,t) and AQ{r,t) = Q'{r,t) - Q{r,t) with Q{r,t) = 
{J-{r,t)) + 'm{Q{r,t)) , and Q'{r,t) the same quantity but in the primed system. 

We next prove that (j71l) admits only one solution, i.e., AA(r, t) is completely 
determined by the averaged dynamics of the current and particle densities, once the 
coupling with the environment via Sen is assigned. To this end we expand (j7ip in 
series about t — Q and obtain an equation for the Z-th derivative of the vector potential 
Ayl(r, t). We thus arrive at the equation 



no(r)(Z + l)AA,+i(r) 



Y^{k + \)ni-k{r)AAk+i{r) 

AQi{r) 
I 



+ ^/i-fe('^) X [V X AAkir^ (72) 

where, given an arbitrary function of time /(r, t), we have defined the series expansion 
_ 1 d'f{r,t) 



Mr) 



(73) 



To prove that the right hand side of ([7^ does not contain any term AAi+i(r) we use 
the fact that the dynamics of any ensemble-averaged operator is given by (j38p . Indeed, 
this implies that the l-th time derivative of any operator can be expressed in terms of 
its derivatives of order k < I, time derivatives of the Hamiltonian of order k < I, and 
powers of the operators Sa and S*^ . The time derivatives of the Hamiltonian do contain 
time derivatives of the vector potential A{r, t), but always of order k < I. Then on the 
right-hand side of ([7^ no time derivative of order / + 1 appears. Equation ([7^ can 
be thus viewed as a recursive relation for the time derivatives of the vector potential 
lS.A{r, t). To complete the recursion procedure we only need to assign the initial value 
AAo(r) — j4'(r, i = 0) — Aexti^, t ^ 0). Since in the unprimed and primed systems the 
densities and current densities are, by hypothesis, equal, the initial condition is simply 
given by nir,t ^ 0)AAo{r) = (*(0)|ip(r, t = 0)|*(0)) - (0)|jp(r, t = 0)|r(0)), 
where jp{r) = (l/2m) J^ilPi^ ^(^ ~ ^0} is the paramagnetic current density operator. 

The same considerations as in [72] about the finiteness of the convergence 
radius of the time series ((7^ apply to our case as well. We rule out the case 
of a vanishing convergence radius by observing that it seems implausible that the 
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smooth (in the ensemble-averaged sense) dynamics induced by (I68p can introduce a 
dramatic explosion of the initial derivatives of lS.A{r,t). If this holds, the expansion 
procedure (j72p can be iterated from the convergence radius time onward. We have 
then proved that (17^ completely determines the vector potential AA(r, t) and thus, 
since Af,xt{T,t) is assumed known, it determines A'{r,t) uniquely, up to a gauge 
transformation. 

To finalize our proof, we consider the case in which U = U' and |^(0)) — 
If this holds, AAo(r) = 0. Then the recursion relation admits the unique solution 
AAi(r) = for any I, and at any instant of time t we have A{r, t) ~ A'{r^ t) (still up 
to a gauge transformation). □ 

The theorem states that the one-to-one correspondence is between the averaged 
current density and the vector potential applied to the system. In view of this, the 
KS vector potential is therefore a functional of the averaged current density, i.e., 
AKs{r,t) = i*'[J(r,i), |*(0))](r,t). Obviously, as in the case of any DFT theory, 
nobody knows the form of this functional, therefore the quality of the results we will 
obtain from our stochastic approach depends on the quality of the approximations we 
are able to make for this functionals. At the moment, this is an open area of research. 
The hope is obviously that standard approximations that have proven very useful in 
the past, like for example the adiabatic- local-density approximation 1131 143] and its 
refinements, provides a solid foundation to explore the reach of this approach. 

By looking at the statements of the theorems, we notice that they all establish 
a one-to-one correspondence between two quantities that have the same "spatial 
dimension". With that we mean that the correspondence is either between a scalar 
(the single-particle density) and another scalar (the potential) or between a vector (the 
single-particle current density) and another vector (the vector potential). It is then 
easy to prove that if one tries to establish a correspondence between the current density 
and the scalar potential, is going to fail 74!. From a mathematical point of view this 
appears quite easy to understand: in creating a connection between two different 
spaces (the one of the densities and the one of the potentials, for examples) one has 
to be sure that the dimensions of these spaces are the same. For this reason, one 
should be able to accept that, if the continuity equation is not valid (see discussion in 
section [^75]) ■ then the only physical information we can extract by using the TDCDFT 
theorem is the current density. For this reason, for example, the results of [T3] appear 
lacking a solid foundation. 

3.2. Closing the system 

The proof we have presented for the theorem 4 leaves open a very fascinating 
possibility. What happens if we do assume that the bath operators Sa are not the 
same in the two systems? Interestingly, it can be seen that the proof remains the 
same, i.e., one arrives at the conclusion that whenever the interaction potentials, 
initial conditions, and now the bath operators are the same then the vector potential 
is uniquely identified [M]. A similar result holds for the scalar potential ^80j. One 
therefore arrives at the amazing conclusion that, if one is interested in the dynamics 
of selected observables then we can close the system, i.e., we set Sq, = in the KS 
system, and still get the exact result. The KS system becomes then a closed and 
isolated system, exactly identical to the ones we already know in standard quantum 
mechanics. We therefore do not need to study the evolution of the system under 
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a random influence of tlie environment. Altliougli tliis possibility opens up a new 
way to look at the dynamics of open quantum systems, care needs to be applied 
in understanding its range of applicability. Indeed, one can easily see that the KS 
scalar or vector potentials needed to reproduce the exact dynamics are functional of 
the coupling with the external environment. A priori, very little is known on the 
way this coupling enters in the KS potentials, and it is very likely that any local 
approximation, so successful in DFT, will fail miserably. Moreover, the KS potential 
in this way loses completely its generality, since it will be strongly affected by the 
system under investigation. It is clear, however, that further investigation along these 
very interesting and potentially promising lines is needed. 

4. Applications of the Stochastic Schrodinger Equation 

In this section, we will discuss some applications of the stochastic Schodinger equation 
to investigate the dynamics of interesting systems. Our excursion must be limited, 
due to the variety of approaches that have characterized this field. Our point of view 
will be starting from systems that evolve according to a SSE as we have discussed in 
section 12.11 This cuts out a certain number of interesting, phenomenological, results 
that do not fit easily in this approach. In selecting the material for this review, we 
have indeed preferred to maintain a certain degree of consistency, rather than present a 
scattered amount of results. The interested reader can look into a series of beautifully 
written books where these topics are covered in more details [H [Ml UHl HI] • 

In selecting the topics, we have also maintained a focus on solid state physics to 
match the audience of this journal: we will discuss spin thermal transport, the onset 
of Bose-Einstein condensation in Alkali gases, the general problem of relaxation in the 
presence of a thermal bath at given temperature. This list cannot be exhaustive of the 
field but, we believe, it gives a flavour of the directions of research one could tackle 
within this working framework. 

4-1- Spin Thermal Transport 

A neat application of a SSE is the investigation of the thermal transport in one 
dimensional spin chains. The simplest system we can think of is a chain of N spin 
1/2 atoms kept at fixed positions. A nearest neighborhood interaction is present, and 
we can consider also the presence of a magnetic field. The two ends of the chain 
are connected to two thermal reservoirs at different temperatures: those two spins 
fluctuate due to the stochastic interaction with the reservoir, energy is transferred or 
absorbed from the neighboring sites, and the whole chain reaches a steady state in 
the long time regime. A simple schematic of the system can be seen in flgure|31 The 
Hamiltonian for this system, when the coupling is not present, is written as 



where, i7„ = {a^^(jy^,a^i), (T*j is the spin operator of the site n along the direction 
i = 2;}, and Q is the interaction constant we have assumed constant through 

all the sites. We have embedded the system into a uniform magnetic fleld which we 
assume oriented along the z axis and of strength lo/2. For a N site chain we have that 



^-2 



N-1 




(74) 



(T^ = 1 «) l(g)s' (g) 1 (g) 1 (g) 1, 



(75) 



n 



N-n-1 
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Figure 4. A chain of spins is connected to two energy reservoirs having different 
temperatures via the spins at the beginning and end of the chain. The reservoirs 
transfer or absorb energy from the chain. The spins are interacting with a nearest 
neighborhood interaction of parameter Q. A position-independent magnetic field 
h = u}/2z is present. 



where 1 is the identity in the 2x2 space and is the i-th 2x2 Pauh matrix. It is 
then clear that each operator cr^ is a sparse matrix of dimension 2^ x 2^. 

An effective way to define the effect of the thermal reservoirs is the following. At 
each instance of time, the direction of the first and last spin are randomized: the new 
spin directions are chosen according to a Boltzmann statistical weight which depends 
on the "local temperature" of the last and first spins, i.e., to the temperature of the 
reservoir they are attached to |81^- An equivalent way is to solve the dynamics via 
either a Lindblad equation or, more conveniently, via an appropriate SSE [811 182j . 
The Lindblad equation for the system is written as 



dt 



m 



where 



VRp{t) 



(76) 



(77) 



E 7;c.i(TL) U^mF} - \FlFip{t) - \p{t)FlFi 

k,l=l ^ 
2 

E 7fe,;(Tfl) (Ckmcj - iGlGipit) - ImclGi) .(78) 

k,l = l ^ ^ 



In ([75)) . the operators F and G are the raising and lowering operators for the first and 
last spin, respectively, i.e., 

Fi = a^, F2 = (7(7, Gi = crj^_i, G2 = o-^_i, (79) 

which are build starting from the raising and lowering spin 1/2 2x2 operators in a 
similar fashion as in (j75p 

s+ = s'= + iP, s" = - iP. (80) 

The coupling matrix 7 is given by ^82j 



7(cj) — A/(a;) 



/(/3,w) + l 



(81) 
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where /(/3,cj) = [exp (/Sw) — 1] is the Bose-Einstem distribution function at 
temperature /3 = l/fc^T, is the bath spectral function, and A is the couphng 

strength, ah evaluated at the energy of the magnetic field, w. Once the matrix 7 is 
diagonalized with ai = and 012 as the eigenvalues, one can rewrite the Lindblad 
equation in the equivalent form 



where the operators and are linear combination of the operators F and G. 
Being in the diagonal Lindblad form, we know that we can easily rewrite the dynamics 
induced by this master equation into an equation for a stochastic wavefunction '^{t) 
[5T] . For this problem, going from the master equation to the SSE offers a great 
numerical advantage. The dimension of each of the matrices entering the Lindblad 
equation is indeed 2^ x 2^, which even for a short chain of 10 spins, means more 
than 10^ elements. Although many of those matrices will be sparse, the same is not 
true for the density matrix itself, therefore one cannot rely on using any efficient 
numerical tool for sparse matrix. In contrast, for the SSE the wavefunction is a vector 
of 2^ elements and since all the spin operators are sparse matrices, one can reduce 
the amount of memory and operations by order of magnitudes, thus making a larger 
number of spins in the chain affordable. Even in this case, the number of elements 
in the chain has to be restricted to about 20. Two ways out are possible. On the 
one hand, a reformulation of the problem in terms of Majorana fermions i83J seems 
to allow for longer chains, up to 100 elements, to be tackled at least to investigate the 
steady-state regime. On the other hand, one could reformulate the problem in terms 
of spin-density waves f43j: by reducing the Hilbert space to the low energy waves, one 
could effectively reduce the dimensionality of the problem. 
In figure [5l the magnetic energy of a given site is plotted, 



where _ff„ = is the magnetic energy of site n, tk — kAt is the discrete-time 

grid of the simulation, and the average in (1551) is obtained from a single realization of 
the random noise using the ergodic theorem where we replace the average over many 
realizations of the noise with the time average of a long time realization. Here, T is 
the number of time steps used for the ergodic average after the steady state has been 
reached. In this simulation, T = 10^. With a thermal gradient, there is an energy 
imbalance between the left and right ends of the chain. The system is therefore kept 
out-of-equilibrium by this imbalance. If we set the two temperatures to the same 
value, then a thermal equilibrium is established in the long time regime |81) . 

In figure, [6] the behaviour of the thermal conductivity of the spin chain as a 
function of the length of the chain is reported. Interestingly, the thermal conductivity 
starts to grow for larger system sizes, suggesting that a ballistic thermal transport 
regime might be reached for this chain. 




i H,p{t) 





k=0 



(83) 
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Figure 5. Ergodic average of the magnetic energy of each site, Hn, over a 
long-time single realization of the stochastic noise. The parameters /S^ = 0.41, 
fiji = 1.39, Q = \ = 0.01, and u) = 1 have been chosen. For the realization a 
final time tfi„ = 10^ in units of uj~^ has been used. From |81l . reprinted with 
permission. 




Figure 6. Thermal conductivity of the integrable spin chain as a function of the 
inverse of the chain length, A'^. For large N it seems the conductivity starts to 
grow suggesting the establishment of a ballistic thermal transport regime. System 
parameters: = 0.25, Pji = 0.5, Q = A = 0.01, u) = 1. Simulation parameters: 
At = 1, tfin = 10'' for Af > 9 or tfi„ = 10® for Af < 8. From [81], reprinted with 
permission. 



4-2. Thermal Relaxation dynamics 

In the derivation of the SSE in section [2TT] we assumed the external environment is in 
thermal equilibrium with temperature T = {ksP)^^ , 

~ Trs e-PtiB • ^^^^ 

From basic thermodynamical considerations [HI [17], we expect the open quantum 
system, which is in contact with this equilibrated heat bath, to evolve in the long- 
time to some steady state that coincides with its thermal equilibrium at the same 
temperature T as the heat bath, 

}T P^^^^ = IT- 
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This dynamics towards thermal equilibrimxi is a non-equiUbrimxi process. In order to 
describe realistic open systems and hence introduce the temperature in the description 
of the dynamics, the equation of motion should ensure, in the long time, this relaxation 
behaviour. If one considers the approximations performed in the derivation of the 
SSE, this begs the question of whether the non-Markovian SSE is still able to describe 
thermal relaxation processes on average. 

In order to understand and investigate these processes, one is interested in the 
conditions for thermal relaxation and how they enter the NMSSE ((23)) . 



id\(l){t)) = Hs\(j){t))dt + \-i{t)V\(j){t))dt 

-iX^V f dTe-'"''^VC{T)\(l){t-T))dt. (86) 



^0 

Here, we consider the environment as described by one operator as this case contains 
all the essential characteristics we would like to discuss. The generalization to an 
environment described by many operators requires some additional care. If the same 
temperature is uniform, we do expect the system to reach a thermal equilibrium at the 
given temperature. If temperature is not uniform, we should not expect any thermal 
relaxation: the system is continuously driven out-of-equilibrium by the temperature 
gradients between the different parts of the environment. One can see that the coupling 
of the system with the environment enters the equation of motion only through two 
quantities, namely the bath correlation function C(t), which is connected with the 
noise, and the system's coupling operator V. Whether the system is driven towards 
thermal equilibrium by the coupling to the heat bath can thus only depend on these 
two quantities. However, one naturally expects that the operator V from the reduced 
Hilbert space of the system does not contain information about the environment, like 
for example its temperature. Hence, thermal relaxation processes have to be highly 
dependent on the structure of the bath correlation function C(r) as it is the only 
quantity that describes the coupling from the side of the environment. 

As we have seen before, the SSE describes an ensemble of wave functions evolving 
under the influence of different stochastic processes. Consequently, only on average one 
will be able to state whether thermal equilibrium is reached and thus we will use the 
in section 12.21 derived master equation ([5^ for the discussion of thermal relaxation 
processes. First of all, we are interested in whether or not thermal equilibrium is 
reached and thus it is sufficient to investigate this with the help of the Redfield master 
equation ((33l) in the limit t ^ oo. One expects, indeed, the master equation ((32)) . 
which corresponds to the NMSSE, and Redfield equation (l33l) to have the same steady 
states. In order that the thermal equilibrium state of the system is a stationary state 
of the non-Markovian SSE, the following equation has to be satisfied 

hm &^ = 0. (87) 

t^oo dt ^ ' 

This requirement and the fact that the equilibrium density operator commutes with the 
system Hamiltonian lead to a condition for the NMSSE to ensure thermal relaxation, 

= ^hm ^ = K(r^V + Vp'^^K^ - VKp''^ ~ (f^K^V + ©(A^), (88) 

where 

/•oo 

k = \^ dTC{T)e-'"^''Ve'^'\ (89) 



CONTENTS 32 

From this one can obtain the conditions for stochastic relaxation processes, to this 
end, we wih change to the energy basis of the system, 

Hs\n) = e„|n),y = ^i;„™|7i)(m|, (90) 
in which ([H5)l can be written as 

7 '■cq 

lim = V \n){l\vnmv^i (91) 
t^oo at ^ — ^ 

X ^ dr|c(T)e-'('^"-^"')^e-'^'^'" + C*(T)e~'(''""'^')^e-'''=- 

- C(T)e-*(''"-<^')^e-''^' - C*(T)e-*('"~^'")^e-'^'4 + 0{\^). 



As discussed before we are interested in system independent conditions under which 
the right hand side of ((9T|) vanishes. Hence, these have to be connected with the bath 
correlation function, 

C(r)=TrB[/5^5^B(r)B(0)], (92) 
and by using the fact that _B is a hermitian operator, one can conclude that 

C*{t)^C{-t). (93) 
As a result of this, the 'half Fourier transform' in ([?T|) can be written as 



+ 1- [ dT(c(T)e-*'^^ -C*(r)e'"^ 
2« Jo I 

= \c{uj)+iD{u:), (94) 

where D{uj) is the imaginary part of this half Fourier transform and as a consequence 
the Fourier transform of the bath correlation function, C(w), is a real- valued function. 
Furthermore, if the half Fourier transform (I94p is analytic in the upper complex half 
plane of io and vanishes faster than as w goes to infinity, one can apply the 

Kramers-Kronig relation, 

D{u) = —P r da^^. (95) 

Here, P J denotes the Cauchy principal value of the integral. 

In the same spirit the half Fourier transform of the conjugate bath correlation 
function can be written as 

1 o r .A--) 



/ dTC*{T)e-"^-' = d{-u) + i—P / da 
Jo 27r 



w — a 



\^{~Lu)+^F{L,). (96) 
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Inserting equations (IMl) and into (lOTj) leads to 



lim — = y \n){l\VnmVml{ {-C{l^nrn) ^ iD{iOnrn)]e^^'''" 



n.rn 




+ (^5(-a;„;) + iF{uj,ni) 

- (ic(-c<.„™) + zFK^)) e-^^" I + 0(A4), (97) 

where Wjj = — ej. We want to point out that in order to satisfy the requirement 
of thermal relaxation, the right-hand side of this equation has to vanish. In addition, 
the Fourier transform of the bath correlation function can be interpreted as the power 
spectrum of the noise and hence describes the probabilities for energy transitions in the 
system. By assuming that this power spectrum satisfies a so-called detailed-balance 
relation, 

C'i^uj) = e^"C(tj), (98) 
(|97p simplifies to 



lim = iX^ y \n){l\snms,ni\ D{ujn,n)e-^'"' +F{< 
t->oo dt ^-^ 



- Diu^mde-f"'' - F{wnm)e-^'- 1 + 0{X^). (99) 

The detailed-balance relation (|98)) ensures that the energy transitions in the system 
are balanced according to Boltzmann statistics. Furthermore, we want to point out 
that (j99p has only imaginary components and by inserting the explicit integrals into 
([Ml) we arrive at 

r dpf iX-^ ^a(a)(l-e-'3(^"-"-))e- 
llm — -f- = — y \n){l\VnmVmi / da< 



t->oo dt 27r ^-^ 



C'(a)(l-e-^('^™-^'-'^)) 



Wmi — a 



OiX^). (100) 



In this expression there is no need to write the principal value anymore since the 
integral is no longer singular. It can be shown that the diagonal components of (jlOOp 
cancel each other, i.e., 

(/|^IO = o. (101) 

This can be done by changing Umn to —Umm changing the integration variable in the 
second integral of (|100p from a to —a, and applying the detailed-balance relation (IMl) 
again. 

As a result, one can conclude that the NMSSE ([25)1 has a stationary solution which 
coincides with the thermal-equilibrium state up to first order in A. Furthermore, if the 
detailed-balance relation is satisfied, the corresponding master equation drives the 
system towards a stationary state that coincides in the diagonal elements in the energy 
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basis with the thermal-equihbrium state up to fourth order. Additionally, neglecting 
either the off-diagonal components of the density matrix in the long-time behaviour 
or the imaginary contribution of the half Fourier transform, 



Im 



dTC{T)e- 



D{lu) w 0, (102) 



the system will be driven towards equilibrium if the equation is considered up to 
fourth order. It can be argued that this imaginary part can be included in the system 
Hamiltonian [12j . leading to the so-called Lamb shift. This means that the equilibrium 
density operator will not commute with this effective Hamiltonian. Nevertheless, this 
energy shift will not introduce dissipative dynamics in the system |84j . 

In conclusion, the bath correlation function should satisfy the detailed-balance 
relation in order to describe thermal-relaxation behaviour. Hence, by the Markovian 
approximation, Ca^i3{T) ~ Da.pS^r), one neglects this property of the bath correlation 
function. Therefore, the description of thermal relaxation processes with the help of 
the Markovian SSE seems to be questionable. 

As an example to demonstrate the previous considerations about thermal 
relaxation processes we will discuss a simple model. Namely, a tight-binding model 
with three sites, 

Hs = -T{c\c2 + c\ci + cjc3 + c\c2), (103) 

where the operators c\ create an electron at the i-th tight-binding position and F is 
the hopping integral. We will couple this electronic system to the electromagnetic 
field inside a cavity. For this model one can derive the power spectrum in dipole 
approximation (c— 1), 

Ceav^tyH = ^ [/(/?, |^|) + 0{-Uo)] if \u:\ < LJ,, (104) 

where /(/3, uj) = l/(exp(/3a;) — 1), 9{uj) the Heaviside function, ujc is a cutoff frequency 
(determined by the dimensions of the system to satisfy the dipole approximation) 
and n is the volume of the cavity. For |a;| > uJc this function vanishes. We want to 
point out that this power spectrum contains naturally the detailed-balance relation 
and hence the model is suitable to describe thermal relaxation processes. 

In addition, one can derive the system's coupling operator for the considered 
model as 

V ^^qj2^-mp)4^p^ (105) 

i,p 

where q is the electron charge and the operator ej creates an electron in the l-th 
energy eigenstate of the system. Here, we have assumed for simplicity that each mode 
of the cavity has the same polarisation direction u, parallel to the chain of the three 
tight-binding sites. We want to point out that the form of this operator should be 
immaterial for the establishment of a thermal equilibrium, since this is only determined 
by the power spectrum p04p . Moreover, if the same operator is used in a Markov 
approximation, a steady state is established that is not the equilibrium configuration. 

In figure [7] we report the probability of occupying the three eigenstates of the 
Hamiltonian (|103p as a function of time obtained from the NMME ([5^ . For this 
dynamics, we have used the parameters, /3 = 1, ajc = l,r = l and A = 0.1. We have 
used a simple Euler algorithm [351 with a time step 5t — 0.001 to numerically solve the 
NMME (1321) ■ Unfortunately, the solution of the NMSSE and master equation requires 
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Figure 7. The probability of occupying the three eigenstates of the Hamiltonian 
II103II as a function of time obtained from the numerical solution of II32II . The 
system is driven towards thermal equilibrium, represented by the three final points 
in the figure, in the long time regime. 

the evaluation of the history integral which contains the correlation function. This 
integral needs to be evaluated at each time step, therefore making the solution of the 
NMSSE unattainable in a reasonable amount of time. This is because one needs to 
build some statistics before performing the averages: for obtaining a reasonable colored 
noise following the discussion in 12 .41 we need to average over thousands of runs. As we 
discuss elsewhere [58 , one can approximate the NMSSE with an equivalent equation 
up to the fourth order in A. This approximation reduces the numerical needs to solve 
the NMSSE, since we arrive at a NMSSE that is local in time. We are able to show 
that, this new SSE (still non-Markovian) and the corresponding master equation give 
the same dynamics. This bring the numerical cost of solving the NMSSE to the same 
level with the Markovian SSE [58] . 

4-3. Ionic motion 

An interesting extension of the result of the preceding section [3] on DFT has been 
proposed recently by H. Appel and M. Di Ventra [S5^. They used the electrons present 
in a certain molecule as the "link" between the ions and the external environment, 
therefore investigating how some vibrational modes, excited at the initial time and 
coupled to the electron motion can be effectively dissipated [S5]. The starting point 
is a theory similar to the one we have discussed in section |3] where one considers 
the total current and charge densities, namely, including in those expressions also 
the ionic contribution. In a system of electrons and ions, where the former are 
connected to external baths, one proves that the total electrical current density is 
in a one-to-one correspondence to the external vector potential, thus establishing a 
stochastic TD-DFT scheme for the system under investigation. This clearly allows for 
energy dissipation from the ionic degrees of freedom via the electrons. In figure |S1 the 
dynamics of the ionic positions of a Neon molecule is shown In the left pane, 

the electrons are disconnected from the environment, and the system undergoes a 
periodic motion determined by the initial conditions. In the right pane, the electrons 
are connected with an external bath that continuously extracts energy from them: the 
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ions relax to the equilibrium configuration of the Neon dimer, with the electrons in 
the ground state. It is clear that this is a simple example of what is an interesting 
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Figure 8. The (classical) dynamics of the ionic positions of a Neon dimer. On 
the left, the system is closed and the dynamics periodic, on the right, the ions 
dissipate energy to the electrons in contact with an external environment. The 
system therefore reaches the ground state after some time. Prom |86| . reprinted 
with permission. 



and promising technique to investigate the correlated dynamics of electrons and ions 
coupled to an external environment. For example, this dissipative dynamics could 
possibly explain the appearance of broadening on the electron and phonon spectrum. 



Bose-Einstein Condensation 

As a last example of application of the stochastic Schrodinger equation we consider 
again the Bose-Einstein condensation [T71 [TS] of Alkali gases [HI]- In this case, we 
depart from the SSE we have discussed in previous sections, and introduce the case 
of particle exchange between the system and its environment. This allows us to 
investigate the dynamics of the condensate formation and predict the condensation 
temperature. To do so, we will need to introduce a stochastic term that is not 
"multiplicative", i.e., it does not multiply the wavefunction as happens in 12.11 but 
"additive" as we will discuss in a moment. 

In the experimental realization of the condensate, the gas is trapped into a 
confining harmonic potential generated by a magnetic field obtained by focusing two 
laser fields [5^1 [SO] ■ In these systems the final temperature is about a few hundreds 
nano-Kelvin. To reach this fantastic result, a novel technique was developed, the 
evaporative cooling of the gas. At temperatures above the condensation temperature, 
the gas follows a Boltzmann statistics. The evaporative cooling selectively removes 
from the gas the particles with the highest kinetic energy. This is possible since 
the highest kinetic energy states are those who have a finite probability density at the 
border of the magnetic trap. Therefore, another laser focused around the boundaries of 
the magneto-optical trap can allow for those bosons to be released from the trap. The 
gas, where the tail of the Boltzmann distribution has been cut, relaxes via boson-boson 
collisions to a new distribution with smaller temperature, since the total kinetic energy 
is lower than before. By selectively removing the hottest atoms, the temperature of 
the Alkali gas can be reduced below the condensation temperature and one can observe 
a condensation both in momentum and position space [SI] [SHI HO] • 

As we have rapidly discussed in section [^75l the Gross-Pitaevskii (GP) equation 
is the standard theoretical tool to describe the physics of the condensate system [ST]. 
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The equation has proven successful in describing both for the dynamics and static 
properties of the condensate. Due to its simphcity and striking successes, many people 
have tried to generalize the GP equation to study the dynamics of the condensate 
formation. In this system, the ground state and the excited states are well separated 
in energy -this separation is dictated by the confining potential- and one could regard 
the excited states as a reservoir of particles and energy for the condensate. With an 
expansion on the particle-particle interaction, by using the Keldish formalism Stoof 
[87] arrived at the following stochastic GP equation 



2 



+ Ve.t{r) - /i - iR{r, t) + g\^[r,tt 

2m 



+ V{r,t). (106) 

Here, describes the condensate phase, Vext is the confining potential, /i is the 

chemical potential, g the interaction constant. Finally, R{r,t) and r]{r,t) are the 
terms coming from the coupling of the condensate to the non-condensate atoms which 
form the "environment" for the condensate. Notice that r]{r, t) does not multiply the 
wavefunction $(r, t), but is added to the equation of motion. For this reason, r] is an 
additive noise for our problem. This is necessary in order for (I106P to allow particle 
exchange with the environment, i.e, the norm of the function ^(r, t) can vary with time. 
Clearly, the term R(r^ t) represents a dissipative contribution denoting the coupling 
between the condensate and the non-condensate phase, while r]{r,t) accounts for the 
incoherent collisions within the gas, analogously to the description of the Brownian 
motion. Interestingly, one can arrive at a similar equation with a completely different 
technique. This technique, developed in parallel and independently from the one 
we present here, looks promising in building a theory to investigate the dynamical 
formation of the Bose-Einstein condensate in Alkali gases [551 Ell EO]- Its first 
application to the investigation of vortex formation close to the transition temperature 
has provided results in good quantitative agreement with the experiments [91]. Here, 
however, we do not feel necessary to expose this theory. The interested reader could 
find a clear exposition in [55] . 

The noise j] is characterized by a Gaussian correlation function 



ij*{r',t')ij{r,t) = -^\r,t)5{t-t')5{r-r'), (107) 

where S'^(r, t) is the Keldish self-energy of the gas of Alkali atoms |87[ 1^155] . As one 
could expect, there is an equivalent form of the dissipation- relaxation relation that 
connects the dissipation term R to the noise 77: 

- 2(lf2/(t)) ' 

where i?(r, e) and YJ^ (r, e) are the Fourier transform of -R(r, t) and YJ^ (r, t) , respectively 
and /(/3, e) is the Bose-Einstein distribution at temperature T. 

The dynamics of the condensate formation can be obtained numerically in one 
dimension. In 2D and 3D, while in principle the equations are still valid, their 
numerical integration is more difficult. Besides the obvious numerical overhead 
associated with dimensionality, this is due to the more complex spectrum of the 
system. Although strictly speaking in one dimension the Bose-Einstein condensation 
for an ideal gas does not exist [TT] [18], in practice one considers the formation of 
the so-called quasi-condensate which corresponds to a macroscopic occupation of the 
ground state of the confining potential. To obtain an effective equation of motion. 
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one splits the wavef unction into a radial and transverse component. Then the radial 
component is assumed to be Gaussian and integrated out. The process gives rise to 
an effective one-dimensional interaction constant, gid- This approximation works in 
an excellent way for the cigar-shaped potentials ^STj where the confinement along the 
radial component is stronger than along the transverse direction. In figurelHlwe report 
the dynamics of formation of the quasi-condensate at different times. At the end of 
the process, a quasi-condensate of 20000 particles is formed at a temperature of about 
100 nK. A key quantity to establish the presence of a phase transition is the so-called 
first-order correlation function, defined as 



5« (-/2, ./2; t) ^ . (109) 

|ci>*(-z/2,i)P mz/2,tW 

This quantity can be measured experimentally [941 1951 I96j . thus providing a test 
for the theory. In figure [9j distances are scaled to the quasi-condensate dimension 




Figure 9. (Upper panel) Snapshots of the particle density during the growth 
process of a one-dimensional quasi-condensate of ^^Na. For the simulations, the 
parameters fi = 30 X Ease where Ease = hui^ with oJz = 27r X 30 Hz for a final 
quasi-condensate ol N = 20000 particles and a final temperature of T = 100 nK. 
Displayed snapshots are at t/teq = 0.05, 0.15, 0.30, and 0.45, where t^q is an 
approximate time for the system to reach dynamical equilibrium. The effective 
interaction g^^ is obtained after averaging over a transverse Gaussian profile of 
width I = ^Jh/ muix where ojj' = 2it X 120 Hz. (Lower panel) The correlation 
function for the quasi-condensate. When the quasi-condensate is forming, leftmost 
plot, the correlation function is expected to decay exponentially with the distance. 
When the quasi-condensate is formed, rightmost plot, the correlation function 
decays as a power law. Reprinted from 93 with permission. 



at equilibrium. The effective radius Rtf{T) varies with temperature starting from 
the Thomas- Fermi expression at zero temperature i?TF(0) — \/2^/nvjj1 if Wz is the 
frequency of the confining harmonic potential in the z direction, and ^ is the chemical 
potential, i.e., the energy of the condensate. For the discussion on how to calculate 
Rtf at finite temperature see ([13]) and references therein. 

The results of figure IHl can be compared with the experimental results obtained 
since 1995, the year of the first experimental observations of the Bose-Einstein 
condensate in ultra-cold Alkali atomic gases [59l [60]. The theoretical modelling is 
able to describe the quasi-condensate formation with a throughout investigation of 
the dynamics. It is also important that the theory contains essentially no fitting 
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parameters, since all of them can be quantitatively determined from the knowledge of 
the experimental setup. 

The Bose-Einstein condensation in Alkali gases has proven to be an important 
tool for a deeper understanding of quantum mechanical phenomena. In particular, 
the relative easiness in manipulating it and its relative stability have allowed for a 
series of interesting experiments on quantum measurement theory to be carried out. 
Indeed, the condensate can allow for either invasive or non-destructive probe of certain 
physical quantities. In this respect, it has been recently possible to develop a stochastic 
theory to describe the case in which the condensate is interacting with a laser field 
which is under continuous measurement. As is the case with the theory of continuous 
quantum measurement in other systems, this theory do not fall into the form of the 
SSE we have extensively presented in this review. The interested reader could start 
from [8l[9]. 

5. Summary 

The stochastic Schrodinger equation is a ubiquitous tool that can be used to describe 
the dynamics of a open quantum system, namely a system in interaction with an 
external environment. The accurate use of this tool requires a deep understanding 
of the physical approximations used to derive, from the Schrodinger equation of the 
total system, an effective equation of motion for the subsystem of interest. In this 
review we have mainly focused only on one class of SSE, those that have a quite 
specific form and structure, derived from a few assumptions on the dynamics of the 
environment and on the coupling between the system and its environment. To give 
to the reader a feeling of how wide is the range of applicability of these concepts, we 
have briefly discussed in section 14.41 a SSE with a different form, in particular with 
an additive noise in contrast to the multiplicative noise we have discussed in the rest 
of this review. Other stochastic equations have been suggested in the past [TU] to 
describe an effective dynamics of the system, without a clear step-by-step derivation 
of the equation of motion. We have not considered that class here. The reason for 
that is simple: each of these models requires a proper introduction, justification, and 
derivation, while we are more interested in presenting a unified, although restricted, 
point of view. The interested reader can start from TD^ and explore the vast literature 
present there. To make connection with the knowledge of the reader, we have clearly 
established the link between the SSE and a master equation for the reduced density 
matrix. From this point of view, one can see the SSE as a "Langevin" type of equation 
and the master equation as its "Fokker-Planck" counterpart. We have shown however, 
that the two theories can lead to significantly different results for physical quantities. 

We have applied the SSE to different systems, from spin thermal transport, to 
Bose-Einstein condensation to thermal equilibrium in electronic systems. The list of 
examples could go further but we will find ourselves in repeating the same, useful, 
concepts. 

Finally, we have discussed an interesting new topic in SSE, i.e., the inclusion of 
particle-particle interactions and its application to describe the dynamics of electronic 
devices. We have shown, in parallel with the DFT theory, how one can build an 
effective description of the open system that in the future should allow us to construct 
an ab-initio modelling of electronic thermal transport, without relying too much on 
any assumption about the properties of the steady-state. 
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The field is still in its infancy. For example, a clear connection between the non- 
Markovian SSE and its master equation counterpart is missing: the actual derivation 
by Gaspard and Nagaoka [27j is incomplete and leaves many questions open, first of all 
the non-uniqueness of the derived equation. A deeper understanding of the condition 
for the establishment of a steady-state or "thermal equilibrium" is also missing. The 
"detailed balance" is indeed not enough to ensure the system is driven towards a 
"thermal equilibrium" , at least in the form known from standard statistical mechanics. 
Also, a thorough exploration of the effect of particle-particle interaction, even at the 
basic level of mean-field theories, is lacking. We therefore foresee an interesting future 
for these techniques. 
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